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We calculate the baryon chemical potential (/is) dependence of thermodynamic observables, i.e., 
pressure, baryon number density and susceptibility by lattice QCD using the canonical approach. We 
compare the results with those by the multi parameter reweighting (MPR) method; Both methods 
give very consistent values in the regions where errors of the MPR are under control. The canon¬ 
ical method gives reliable results over = 3, with T being temperature. Multiple precision 

operations play an important roll in the evaluation of canonical partition functions. 
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I. INTRODUCTION 

Quantum Chromodynamics (QCD) is the fundamental 
theory describing the strong interaction. It is well known 
that QCD has the rich phase structure at finite tempera¬ 
ture and density [T]. And yet, regions that we can access 
with the perturbation are limited. Currently, the most 
promising method to explore the phase diagram is lat¬ 
tice QCD simulation which is first principle calculation 
of QCD. 

Although the lattice QCD simulations are very success¬ 
ful to analyze the phase diagram of a finite temperature 
system, at finite density they have a severe problem, so- 
called sign problem, and outcomes of the first-principle 
calculation can be available only at small chemical po¬ 
tential range. In finite temperature and density systems, 
lots of physically interesting targets such as the early 
universe, neutron stars, quark matters are waited to be 
explored. Therefore, it is quite desirable to explore meth¬ 
ods for investigating finite density QCD systems from ab 
initio calculation; this is the one of the urgent subjects 
in particle physics and nuclear physics. 

The canonical approach we study in this paper is a 
promising candidate for this purpose. In Ref. [2], the fu- 
gacity expansion by a method of the hopping parameter 
expansion was constructed as a winding number expan¬ 
sion, and the chiral condensate as well as the thermody¬ 
namic quantities is studied. More detailed analyses were 
performed in Ref. [3] in a wide range of the temperature 
and chemical potential regions, and an indication of the 
transition was first observed below Tc and finite baryon 
density. In this paper, we address two questions: 

1. Does the lattice canonical approach produce con¬ 
sistent results with the MPR ? 

2. In obtaining the canonical partition functions for 
large baryon number, what is a role of the multi 
precision calculations ? 


Basic concept of the canonical approach in QCD 

In Nf flavor QCD case with the degenerate quark 
masses, the grand canonical partition function at finite 
temperature T and finite quark chemical potential fiq is 
given in the path integral formalism as follows. 

Zgc{T,^i,) = J d[U]{detA{^i,)}^fe-^^, ( 1 ) 

where detA(/ig) is the one flavor fermion determinant 
and Sg is the gauge action. Because the fermion deter¬ 
minant has the property 

[det A(/iq)]* = detA{-n*g), (2) 

the Monte Carlo measure {det becomes 

complex number at finite real chemical potential and 
the standard Monte Carlo method breaks down. Conse¬ 
quently, we cannot study finite density thermodynamics 
with standard grand canonical method. This difficulty is 
called sign problem. 

A system described by the grand canonical partition 
function Zgc{T, jig) is equivalent to a system described 
by the canonical partition function Zc{n,T) with fugac- 
ity in thermodynamic limit. The relation of two 

ensembles can be written as a fugacity expansion [4] us¬ 
ing eigen vectors of number operator N \n) = n\n), 

ZGc{T,fiq)=Tre-^^-f^'>^y^ 

OO 

= Y, («l |n) 

n= —OO 
OO 

= ^ Zc{n,T)e^>^^/^, (3) 

n= —OO 

where is fugacity. If we have the canonical parti¬ 

tion functions, Zc{n, T), for all net quark numbers n, we 
can construct the grand canonical partition function as 
a polynomial of fugacity with coefficients Zc- From this 
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formula, one can obtain Lee-Yang zeros [5], which reflect 
the system’s critical nature [6]. 

The canonical partition functions are constructed 
through the Fourier transformation of grand canon¬ 
ical partition function at pure imaginary chemical 
potential [7], 


Here n,m are space-time coordinates on a lattice, n is 
hopping parameter and fiq is the quark chemical potential 
which is introduced to the temporal part of link variables. 

In order to obtain canonical partition functions, we 
need to compute the grand canonical partition functions 
at various pure imaginary chemical potential values in 
Fourier transformation Eq. 

We use the reweighting method to evaluate the grand 
canonical partition function. 


where fii G M. Eq. (© tells us that the fermion determi¬ 
nant is real in the case of pure imaginary chemical poten¬ 
tial. Monte Carlo simulations can then be performed and 
the canonical partition functions are obtained by Eq.Q. 
Eq. also insists that the canonical partition factions 
are real number because the grand canonical partition 
function is even function (charge conjugation invariant) 
in terms of chemical potential. Considering this feature 
with Eq.(|^, one can find that canonical partition func¬ 
tions are real and positive also in the context of canonical 
approach. 

Once Zc are available, we can construct the grand 
partition function by Eq. (§ at any real quark chemical 
potential. This is because the chemical potential depen¬ 
dence of the grand canonical partition function appears 
only through fugacity, which is the variable of the 

polynomial, and not in the coefficients, Zc in Eq. (§ ; the 
effect of the chemical potential appears through the fu¬ 
gacity and the canonical partition function plays just a 
role of coefficients in the fugacity expansion of the grand 
canonical partition function. 


II. FRAMEWORK 

A. Winding number expansion of grand Partition 
function 


In this work, we employ the RG-improved gauge action 

/3 




CO E E 


n,iJ,<u 


n,iJ,<u 


( 5 ) 


with Cl = —0.331 and cq = 1 — 8ci, and the clover im¬ 
proved Wilson fermion action with the quark matrix 

A(77/, 771, /7g) =Snm SW^nm ^ ^ 
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+ (1 + 7i) U} 

-/c[e+'^^“(l-74)!74(n),5„,„+4 


— 1 


(6) 




det A(i/i/) 
det A(/io) 

detA(i/i/) 
det A(/io) 


1^/ 


Mo 


^gc(mo), 




( 7 ) 


where /xq = 0 or pure imaginary values. We can then 
evaluate the canonical partition function as 


Zc{n,T) 

Zc{0,T) 




(8) 


We adopt, here, normalized canonical partition functions 
Eq.(8) in order to avoid the extra constant Zccif^o) in 
Eq.(7); this step does not affect the physical result. Now 
the evaluation of the grand canonical partition function 
is reduced to the calculation of ratios of fermion deter¬ 
minants in Eq.Q. 

Performing the hopping parameter expansion in the 
logarithm, we write the fermion determinant as 


det A(i/i/) = exp Trlogjl — /^A(i/i/)} 


= exp 


-TrE-Q7w) 


i=i 


( 9 ) 


The trace is taken over space-time, spinor and color. 
Here, we used the following identity for arbitrary matrix 


A: 


= = ( 10 ) 

and log is expanded assuming k is small (hopping param¬ 
eter expansion). 

The contribution of the trace in Eq. § comes from 
all closed loops on a lattice, and the chemical poten¬ 
tial dependence comes from specific closed loops winding 
along positive and negative time directions. We can thus 
classify the trace of the hopping parameter expansion 
in Eq.(|^ according to the winding number which is the 
number of net windings along the time direction. As a 
result, we can reach following expression with coefficients 
Wn and complex fugacity Here n represents the 

winding number. 


det = exp 


E 


Wr,e 


infj,! jT 


( 11 ) 


We call this expression as ‘winding number expansion’. 
The negative winding number appeared in Eq. 0 stands 
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for the winding along negative time direction. The co¬ 
efficients Wn has no chemical potential dependence; the 
chemical potential dependence appears in the fugacity. 
Consequently, we have only to calculate Wn from given 
gauge configurations to obtain grand canonical partition 
functions at desired pure imaginary chemical potential. 


B. Constraint on canonical partition function from 
symmetry of QCD 

Roberge and Weiss pointed out that the QCD grand 
canonical partition at pure imaginary chemical potential 
has the following periodicity [8] 


where Ng = = Ny = and T~^ = Nfa with a lattice 

spacing a. The deviation of the pressure from = 0 is 
given by 


Ap{fiB,T) 


p(0,T) 


J 14 


J14 



3 

log 


V Zgc{0,T) )' 


(18) 


The dimensionless baryon number density tib/T^ and 
susceptibility x/T^ are 


UBi/J^B^T) _ d p{pb,T) 
T3 - d{pB/T) T4 


(19) 




ZgG = ^GG 


ipi 


2'Kik\ 

^ J 


( 12 ) 


where /c G N. Using Eq.( 12), we rewrite the grand canon¬ 


ical partition function as 
2 


p/i/ \ _ 1 ! ipi 

GC y 2 ^ j ^ / V GG ( rp 


1 


/e=0 


T 


) 


(13) 


x{^b,T) p{^b,T) 

T2 d{i2B/Ty T4 


( 20 ) 


III. NUMERICAL RESULTS 
A. Lattice set up 


Then, we get the following relation. 


Zc{n,T) 



We obtain the following important constraint on the 
canonical partition functions, 

ZG{ny^3k)=0. (15) 

Note that this holds both in the confinement and the 
deconfinement phases. 

Now the grand partition function can be written as 

00 

Zgc(T,mb)= E (16) 

B= — oo 


where B e N. Because this quantum number B can be 
interpreted as net baryon number, pB can be regarded 
as baryon chemical potential which is related to quark 
chemical potential as /i^ = 3pq. 


C. Thermodynamic observables 


In a homogeneous system, the dimensionless equation 
of state at (pb^T) is given by 


p{Pb,T) 

J14 


1 


WT3 


log (mb, ^) 


= ( ^ ) logZGG(MB,T), (17) 


We adopt 2-fiavor clover improved Wilson fermion ac¬ 
tion with Csw = (1 ~ 0.8412//3)“^/^ evaluated by one- 
loop perturbation theory and RG-improved gauge action. 
All simulations were performed on x NyXN^xNt = 
8 X 8 X 8 X 4 lattice. We considered /3 = 2.00, 1.95, 1.90, 
1.85, 1.80, 1.70 which correspond to T/Tc = 1.35(7), 
1.20(6), 1.08(5), 0.99(5), 0.93(5), 0.84(4)0. The val¬ 
ues of hopping parameter k was determined for each P 
by following the line of the consistent physics in case of 
vriT^lvrip = 0.8 in Ref.0. 

We generated gauge configurations at /io = 0 with the 
hybrid Monte Carlo (HMC) method. The step size dr 
and number of steps W of HMC were set to dr = 0.02, 
Nn = 50 so that the simulation time was dr x = 1. 
After the first 2000 trajectories for thermalization, we 
adopted 400 configurations every 200 trajectories for each 
parameter set. 


B. Instability of Fourier transformation in 
canonical approach and its solution 

Before proceeding to our numerical results, we refer 
to the instability of Fourier transformation in canonical 
approach and then discuss our strategy to avoid it in this 
subsection. 

Since the fugacity expansion of grand canonical parti¬ 
tion function, Eq.(|^, converges at real baryon chemical 
potential, the canonical partition function Zn must be¬ 
come smaller when the net baryon number |n| becomes 
large. This means that we have to deal with quite small 
values as the results of Fourier transformation. This step 
is quite difficult in the point of view of numerical calcula¬ 
tion because the Fourier transformation is an oscillatory 
integral. 
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FIG. 1. (Color online) Cancelled significant digits in calcu¬ 
lations of DFT at over Tc (upper red points) and below Tc 
(lower green points). 


1. Instability of Fourier transformation 


In order to reduce the effect of this cancellation, we 
should increase significant digits as a solution. Let us 
consider the next calculation with 22 significant digits. 

1.234567444444444444444 - 1.234566111111111111111 
(22 significant digits) 

= 0.000000133333333333333 . . 

(16 significant digits) ^ ^ 

Although six significant digits are definitely lost in this 
calculation, 16 significant digits survive in the final result. 

Summarizing above, the precision of the result can be 
kept by increasing significant digits of variables in this 
way. Figj^ represents the cancelled digits in the calcu¬ 
lation of Zc{B^T) with 16, 32, 48, 64 precision calcula¬ 
tion. According to this figure, for evaluating Zc{B,T) 
to larger n, it is essential to increase significant digits of 
variables. 


In numerical calculation, Fourier transformation 
Eq. is computed by discrete Fourier transformation 
(DFT) as 


C. Thermodynamic observables at finite real 
baryon chemical potential 


Zc{n,T) 


1 

N 


N-l 

k=0 



( 21 ) 


where N is the interval number of DFT. Because DFT 
is just a discretized version of Fourier transformation in 
continuum theory, the instability of DFT in canonical 
approach is simply caused by the numerical errors. They 
are classified into rounding error, truncation error, can¬ 
cellation of significant digits and loss of trailing digits. 
The instability of DFT does not come from truncation 
error since DFT is not infinite series. Accordingly, it is 
quite natural to consider that the instability originates 
from cancellation of significant digits or loss of trailing 
digits or both of them. In this work, we actually moni¬ 
tored the behavior of all variables in our DFT program 
in order to study the effect of these two errors. As the 
result, we found that cancellation of significant digits is 
not negligible in DFT program. Fig{l] represents can¬ 
celled digits in DFT. For example, 80 digits are cancelled 
in case of ^ = 1.80, B = n/3 = 40. We also found that 
the appearance of the cancellation does not depend on 
temperature of a system. 


1. Calculation procedure 

First, we computed coefficients of winding number ex¬ 
pansion Wn up to n = 120 with 400 configurations in 
all temperature cases. We used 64 (above Tc) and 128 
(below Tc) noise vectors to calculate the trace in the 
fermion determinant Eq. Then, we evaluated the 
grand canonical partition functions at various pure imag¬ 
inary chemical potentials using the winding number ex¬ 
pansion with Wn- After that, we evaluated the canon¬ 
ical partition function through Eourier transformation 
and thermodynamic observables. We adopted multiple 
length precision calculation [10] with 400 significant dig¬ 
its in order to keep the sufficient precision except for the 
calculation of gauge configurations and Wn- Gauge con¬ 
figurations and Wn were computed with double precision, 
i.e., around 16 significant digits. 

Note that canonical partition functions are complex 
number in numerical calculations because of numerical 
errors. Therefore, we adopt the only real part of the 
canonical partiton function. If the real part of the canon¬ 
ical partition function is negative at some baryon number 
ub fo the first time, we adopt the result up to — 1 as 
canonical partiton functions. 


2. Solution of instability - multiple length precision 
calculation 


2. Estimation of truncation error in fugacity expansion 


Cancellation of significant digits arises from the follow¬ 
ing types of calculation, 

1.234567- 1.234566 = 0.000001 .^2) 

(7 significant digits) (1 significant digit) ^ ' 

In this case, six significant digits are lost. 


In numerical calculation, we have to deal with the fu¬ 
gacity expansion of the grand canonical partition func¬ 
tion as the finite series 

Zgc{T,^b)= ^ (24) 
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FIG. 2. (Color online) Relation between behavior of Zc(B,T) and precision of variables, at over Tc (left panel) and below Tc 
(right panel). Both panels are plotted with some precision at 16 digits (double precision, uper red points), 32 digits (second 
green points), 48 digits (third blue points), and 64 digits (lowest cyan points). 


Therefore, we have to judge the baryon chemical poten¬ 
tial region where results are free from the truncation er¬ 
ror. There may be several possible ways to analyze the 
effect of the truncation error; In this work, we use the 
following analysis. 

First, we evaluate the expectation values of thermody¬ 
namic observables (0(/iB))7v with Eq.(24). Next, we 
calculate the expectation values 

tracting one from A^max in Eg. ([2^. After that, we evalu¬ 
ate the relative error from these two expectation 

values; and in this work we judge that the expectation 
value is reliable if the relative error is less than 10“^, 


RohifJ^B) = 1 — 


(O) 


^max 1 


(O) 


JV„ 


< 10“^ 


(25) 


In this way, we can ensure that expectation values of ther¬ 
modynamic observables in the baryon chemical potential 
region determined by above analysis have two significant 
digits at least against the truncation error. 


3. Thermodynamic observables 

Using the error estimation described in the previous 
subsection, we analyze the chemical potential dependence 
of thermodynamic observables and study the validity of 
canonical approach. Eirst, we examine the pressure. 
Eig|^ shows that the results of pressure above Tc do not 
suffer from large error up to around /j^bIT = 5, and the 
results below Tc are under control upto/i^/^ = 3.5 — 4. 
On the other hand, the result just below Tc case is reliable 
only up to around {Ib/T = 3. This may be because we 
make configurations at /io = 0 and the configurations are 
suffered from the fluctuation caused by the phase tran¬ 
sition at zero density. We may get clearer signals if we 
generate conflgurations at pure imaginary chemical po¬ 
tential because Tc at pure imaginary chemical potential 
is higher than Tc at zero chemical potential. 

In Eigj^ we see that pressure calculated by canonical 
approach are consistent with results by MPR. 



FIG. 3. (Golor online) Ghemical potential dependence of 
pressure. Red, green, blue, cyan, magenta and brown points 


are the results at T/Tc — 1.35,1.20,1.08,0.99,0.93 and 0.83. 
Upper bound of baryon chemical potential is determined by 
Eq.([25). 
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FIG. 4. (Golor online) Gomparison of pressure calcu¬ 
lated by canonical approach and multi parameter reweight¬ 
ing method. Basically, color of data are same as Figi 
Extra color: dark-red, dark-green, dark-blue, dark-cyan, 
dark-magenta and dark-brown points are the results at 
T/Tc = 1.35,1.20,1.08,0.99,0.93 and 0.83 calculated by 

multi-parameter reweighting method. 
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FIG. 5. (Color online) Chemical potential dependence of 
baryon number density. Color of data are same as Figi 
Upper bound of baryon chemical potential is determined by 
Eq.(I^). 


FIG. 7. (Color online) Chemical potential dependence of 
susceptibility. Color of data are same as Figi Upper bound 
of baryon chemical potential is determined by Eq.(25). 



EIG. 6. (Color online) Comparison of the baryon number 
density calculated by canonical approach and multi parameter 
reweighting method. Color of data are same as Eig|4] 

Next, we consider the expectation value of the baryon 
number density. In Figj^ we find that the results are 
reliable up to around (Ib/T = 4 = 3 — 3.5) above 

Tc (below Tc). While, the reliable baryon chemical po¬ 
tential range of the result just below Tc is limited up to 
IIb/T = 2.4. This may be caused by the same reason in 
the analysis of the pressure. 

Figi tells us that the canonical approach is consis¬ 
tent with MPR method also in the baryon number den¬ 
sity case. Moreover, we observe that the gradient of 
the baryon number density, as a function of baryon 
chemical potential becomes smaller as the temperature 
decreases. In zero temperature case, is expected to 

be zero up to (Ib/T = vtibIT^ where is the lightest 

baryon mass of the system and it becomes a finite value 
at this point. Indeed, the data at T/Tc = 0.84 shows 
such a feature. 

Finally, we investigate the susceptibility. Figi shows 
that the results above is reliable up to around {Ib/T = 
3.5, while the results below Tc is reliable up to Tc = 2.4 — 
2.9. From Figj^ we find that canonical approach is very 
consistent with MPR method also in the susceptibility. 


The susceptibility as a function of the /i^/T does not 
show a clear peak; we do not see yet the signal of the 
phase transition between confined phase and deconfined 
phase. 


IV. SUMMARY AND OUTLOOK 

In this paper, we find that the canonical approach is 
consistent with MPR method. Moreover, the canonical 
approach provides reliable results beyond i^bIT = 3 in al¬ 
most all observables. This is very encouraging for the first 
principle calculation of finite density QCD, because other 
methods such as MPR method, Taylor expansion method 
and imaginary chemical potential method give reliable in¬ 
formation practically only up to /i^/T = 3. The multiple 
precision calculation greatly contributes this conclusion. 

The canonical approach has been investigated several 
times miiiHis] ; We brush up the method here and find 
that it is a useful and promising method. But, we have 
to improve our method further to obtain results in more 
realistic condition, i.e., lighter quark mass, large volume, 
finer lattice spacing and larger density. Although the 
hopping parameter expansion gave very interesting re¬ 
sults as we saw in this paper, the next step is to calculate 
the fermion determinant without this approximation; we 
have learned in this paper that the key point is to cal¬ 
culate the determinant at imaginary chemical potential 
values that are Fourier transformed with high accuracy in 
Eq.Q. This requires more computational resource than 
the work reported here, but within scope of the next gen¬ 
eration high performance era. 
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FIG. 8. (Color online) Comparison of susceptibility calculated by canonical approach and multi parameter reweighting method. 
Color of data are same as FiglU 
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